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ABSTRAC7 

A doubling principle first used by van de Hulst has been de- 
veloped to provide rapid and accurate results for the problem of 
diffuse reflection from a plane parallel atmosphere. The method 
described here eliminates the need of numerically solving an equa- 
tion of radiative transfer by beginning with a layer of such small 
optical thickness (t~2 that the initial scattering and trans- 

mission functions are given by the phase function. An application 
is made to spectral scattering by clouds in the near-infrared where 
the phase function is strongly peaked in the forward direction. 


I . INTRODUCTION 

The problem of radiative transfer in plane parallel atmospheres 
has been solved "exactly" only for very simple phase functions so 
that a numerical approach is required for most practical applica- 
tions. van de Hulst was the first to note that considerable com- 
puting time could be saved if the transfer equation was solved 
only for a layer of small thickness t q and this solution was then 
used to obtain solutions for layers of thickness 2t q , 4 t q , etc. 
by a "doubling" procedure. Numerical verification of the ef- 
ficacy of this method was made by van de Hulst and Grossman (1968) 
who used the Neumann series method (iteration in orders of scat- 
tering) to obtain solutions at t q . Then, after several doublings, 
they fit their numbers to asymptotic (t -*• ») results. 
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In this paper we describe a computing method which gives fast 
and accurate results for all optical thicknesses. We begin with a 
layer of such small optical thickness that the scattering and trans- 
mission functions are given to high accuracy by the phase function 
for single scattering. This accuracy allows the doublings to be 
carried to large optical thicknesses and asymptotic solutions are 
not required. With this method it is not necessary to solve the 
usual transfer equation and the computing program is applicable to 
all phase functions. 


II. DOUBLING EQUATIONS 


In the " standard"problem of diffuse reflection from a plane- 
tary atmosphere considered by Chandrasekhar (1960) a parallel beam 
of radiation of net flux ttF per unit area normal to itself is in- 
cident from the direction ( - u 0 ,<f> 0 ) on .a plane parallel atmosphere 
of optical thickness x. Now if the scattering function S(x;y,*;y o ,<}> o ) 
and the transmission function T (x ; y , <f> ; y Q , <f> Q ) (defined by eq. [l20] , 
p . 20 , Chandrasekhar 1960) are known then S (2x ; y , * ; y Q , <l> 0 ) and 
T(2x;y,4>;y o ,<i> o ) may be found by adding the ways in which diffuse 
radiation may escape from the double layer: 


S(2x;iif4>;Uor$n)F _ S(T;y,4 >;yo,<j>o)F + T / y — — y, (x;u,d>?y ,4> )F e T ^ y o 

4y * 4y u 0 00 


y o 


tt y 


2tt 


T(T,•y,<()?y , E 0 (T;y' / 


• y )dy 'd<t> 'Fe t/m o 

0 o 


( 1 ) 



-4- 


+ e " T/l, 3^r 0r r o (T,u ' * '■ u ’ ' * ’ ) T(T ' u ’'^i‘ i °'*° ) dp -at- F 


1 

f 1 

f 2 ’ 

1 A 

4 ir y 

Jo- 

L J 

O' 


2 TT 


T ( x ; y , 4> ; u " , 4> " ) j^r I 0 (t ; y " , <fr" ; V ’ ,♦ ’ ) 


dy * d 4> 1 dy "d<j> " F 


where 


E 0 (T;y^;y o ,^ o ) = J 

n= i , J 


S n (x;y,4>;y o ,4> 0 ) » 


( 2 ) 


S 1 (x;y,<)>;y o ,<J> o ) = S ( t ; y , <t > ; y Q , $ 0 ) , 


(3) 


and 


S n ( T ; y , <J> ; y Q , <p Q ) 


1 

37 


2 7T 


v du 1 

S ( t ; y , 4> ; y ' , 4> ' ) s n _]_ (x;y' ,<J>' ;y o ,<t» o ) — jp- 


d*'. (4) 


The factor F/4y is included in equation (1) to make clear the 
physical basis for that equation. The first term on the right 
side is the intensity of radiation diffusely scattered from the 
upper layer without interaction with the lower layer. The re- 
maining terms are the radiation transmitted downward by the upper 
layer (diffusely or without interaction) , scattered any number of 
times back and forth between the two layers , and finally transmitted 
upward by the upper layer. Similarly an equation for T (2x j y , <J> ; y Q , <j> Q ) is 
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e 00 n=2 , 4 . . . n 0 


( 6 ) 


The above equations are equivalent to those of van de Hulst (1963) . 
We have further derived equations which include the effect of a re- 
flecting planetary surface and internal sources (e.g., thermal emis- 
sion) but we will omit these since they easily follow from the 
doubling principle. 
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III. INITIAL LAYER FOR COMPUTATIONS 


The above equations provide the scattering and transmission 
functions for layers of increasing thickness provided those func- 
tions are known for any initial layer. However, accurate solutions 
of the radiative transfer equation are often difficult to obtain 
and inaccuracies can lead to an increased error after many doub- 
lings. But for very thin layers the radiation singly scattered is 
an approximate solution which is known for any given phase function. 
For t<<1 the ratio of multiply scattered radiation to singly scat- 
tered radiation is =x and hence the relative error in S or T is 

7 _-25 

only about 1 part in 10 at an initial optical thickness - 2 

Moreover, the above equations suggest that the relative error will 

not increase much during the doublings until t approaches unity 

(since for t<<1, S and T are <<1) and the numerical computations 

verify this. Consequently, if the computations are started at very 

small x , then when t reaches values ~1 the scattering and trans- 

0 

mission functions are sufficiently accurate for the doubling to be 


carried to very thick layers. 

The well known expression for the intensity of singly scat- 
tered radiation (Chandrasekhar 1960) gives us the following initial 

functions 
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T(t o ;w,*;u o ,<|> 0 ) 
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exp(-t o /y o )-exp 


(-T o /p)|p(-y,<j>;-u o f<l , <:) ) (8) 


where P(y,4>?u o ,4> o ) is the phase function normalized to uj o , the 
albedo for single scattering. These equations provide the initial 
values for S and T but since the terms in brackets involve the dif- 
ferences of nearly equal numbers some of the accuracy would be lost 

on a computer; hence it is better to expand the exponentials in 

2 

power series and keep terms through t o - 

The sums (eq. [2] and eq.£6]) may be truncated at some n=N~5 
with the omitted terms replaced by the geometric formula. This is 
a result of the fact that radiation scattered back and forth between 
the two layers tends toward an isotropic distribution after a number 
of such scatterings and hence the ratio of successive terms in the 
infinite series approaches a constant value (a conclusion already 
reached by van de Hulst (1963) for thick layers) . 

Considerable advantage is gained in computing time if the azi- 
muth dependent functions are expanded in Fourier series in 
If the phase function may be expanded in cosines of the scattering 

angle then 


S ( t ; v , <j) ; v Q , <t> Q ) 
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S m ( t ; y , y o ) cos mU-<j> 0 ) 


(9) 
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where 


„m , v 1 
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= 0 if m ? 0 


and 


0 ,m 


= 1 if m = 0 


and similar expressions exist for T. From these it is easily 
verified that each component doubles independently, i.e., each 
pair S m and T™ satisfies equations like (1) and (5) but the 
terms in the infinite series are given by 


_m , 

S, (t ?y ,y ( 


) = S m (t ;y ,y q ) 


( 11 ) 


and 


s ™(T;y,y o ) 
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The computational method is straightforward and efficient 
since it consists of evaluating sums and integrals which may be 
replaced by sums through Gauss quadrature. The same program is 
applicable to any phase function but for elongated functions an 
increased number of Gauss divisions must be used to maintain the 
For finite x the results check exactly to the 6 digit 


accuracy . 
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accuracy available for a few simple phase functions and the same 
is true when the results for t = 2 5 are compared to the "exact" 
values for T = « . In the azimuth dependent case with elongated 
phase functions the results check to the 4 digits of accuracy 
which is estimated for calculations by Hansen (1967) , who used 
the invariant imbedding approach, and for calculations by Potter 
(1968) , who used the Neumann series method. 

Except for strongly peaked phase functions the computing 
time is negligible. For aerosol phase functions varying by 3 
orders of magnitude from their peak value to their low value 
(requiring ~ 50 terms in the cos m<j> expansion) the total computing 
time is 3 minutes on the IBM 360/95 for the scattering and trans- 
mission functions (and derived quantities) for every (not each) 

T multiple of 2 from 2 -25 to 2 7 for 13 values of y , 13 values of y 
and any reasonable number of values of <t> — 4> Q - 

IV. REFLECTIVITY OF CLOUDS 

To illustrate the utility of the computing method the diffuse 
reflection from ice clouds in the near-infrared has been computed. 
The phase functions including the albedo for single scattering, 
u were kindly computed by H. Cheyney from Mie scattering theory 
for spherical particles. The "cloud model" size distribution of 
particles (Deirmendjian, 1964), which has its maximum at diameter 
8y, was used and the optical constants were taken from Irvine and 
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Pollack (1968) . Although natural ice crystals are highly non- 
spherical , it is often assumed that randomly oriented nonspherical 
particles may be approximated by spheres. The results for water 
clouds would be qualitatively the same as those for ice since the 
optical constants of water and ice are similar in the near-infrared; 
however , the corresponding features would be shifted slightly (<.ly) 
in wavelength and in absolute value. 

The calculations were made at 20 wavelengths between .95y and 
3.6y including each wavelength at which the tables of Irvine and 
Pollack show the absorbtivity to be an extremum. The computed 

values of (o (A) are given in Table 1 along with the relative optical 

0 

thickness of the cloud as a function of wavelength. Representative 
phase functions shown in Figure 1 illustrate the increasingly sharp 
forward scattering toward shorter wavelengths and the damping of 
backscattering at wavelengths (~3y) where the absorption is large. 

The spherical albedo as a function of wavelength is shown in 
Figure 2 for several optical thicknesses of the cloud. Although 
the scattering is not far from being conservative (6i Q >.95) for 
A<2.7y, it is obvious from the albedos for increasingly thick 
layers that multiple scattering greatly enhances the "absorption" 
features near 1.5 and 2.0y; the minor features at 1.03 and 2.6y, 
however, are just visible for thick layers. Since dep® nc *s 

strongly on particle size it will be worthwhile to make calcula- 
tions for other particle size distributions and to compare the 
results to measurements of atmospheric and laboratory clouds. 
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Angular distributions of the scattered light are shown in 
Figures 3-8 for the phase function for X=.95y. In the calcula- 
tions for those figures the scattering was assumed to be conserva- 
tive (w =1) ; hence the results are approximately valid for clouds 

0 

0 

of either water or ice spheres for wavelengths from 4000A to 
lGyDOOA because the absorbtivity of both materials is negligible 
in that region and the real refractive indices are within the 
small range 1.30<n r <1.34. The effects of Rayleigh scattering and 
molecular absorption are not included. For wavelengths other 
than . 95y the results apply to a size distribution peaking at a 
diameter d such that d/X * 8/95. In the calculations for Figures 
3-8 140 terms were employed in the cosine expansions and 40 y 
divisions on the interval (0,1); the number of divisions was 
varied and the results suggest that the errors are < 1% for all 
angles. The reflection function which has been graphed is 
defined by 


s (i;y ,<|>;y 0 ,<f> 0 ) 

R(T;y,*;y Q ,* 0 ) = ^ 


(13) 


Reflectivities are shown for thin (t=1) , intermediate (t=4) 
and very thick (x-128) clouds for which cases the spherical 
albedos are 13.3, 34.0 and 93.5%, respectively. Thus the 
diffuse transmission is significant even for extremely dense 
clouds; this is due to the fact that the scattering i» 
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conservative and peaked in the forward direction. 

For conservative scattering it is often assumed that 
thick clouds may be approximated by a Lambert reflector 
(R(t;jj,$;v 0 ,<|> 0 ) « y Q ) in which case the reflectivities would 
be horizontal lines in Figures 3-8. The graphs confirm a 
trend toward Lambertian reflection but Figures 3-5 indicate that 
even for t=128 there is still a nonnegligible dependence of the 
reflectivity on y , especially when the incident direction is 
near grazing (cos _1 y o ~ 90°). Figures 6-8 verify that much of the 
azimuth dependence of the reflectivity disappears for thick lay- 
ers but this also does not hold for large incident angles. 

In Figures 6-8 the intensity of light singly scattered from 

the cloud is also shown; these curves are only given for the 
layer of thickness t= 1 since the increase for thicker layers is 
small. Except for thin layers and near grazing incident and 
emergent angles it is apparent that the light singly scattere 
represents only a small fraction of the total reflected radiation, 
and it is for this reason that the features of the phase function 

are much diluted for thick layers. 

The results of this section suggest that it may be possible 
to obtain a certain amount of information about the phase function 
of cloud particles from measurements of the angular distribution 
of reflected visual light, but since measurements near the graz- 
ing direction are difficult that approach may not be very 
practical except for thin clouds. However, for thick clouds the 
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spectral features in the near infrared become strong and hence 
that wavelength region appears especially promising for cloud 
identification. 
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TABLE 1 

SINGLE PARTICLE ALBEDO AND RELATIVE CLOUD OPTICAL THICKNESS 




.95 

1.03 

1.10 

1.30 

1.52 

1.70 

1.85 

1.90 

2.00 

2.15 



T ( X ) 


u> 


99995 

99986 

99994 

99929 

97126 

99201 

99754 

98574 

94775 

97887 


1.0000 

1.0057 

1.0096 

1.0223 

1.0350 

1.0445 

1.0520 

1.0545 

1.0597 

1.0688 


2.30 

2.50 

2.55 

2.60 

2.70 

2.80 

3.00 

3.20 

3.40 

3.60 


o <»> 

T ( X ) 

99424 

1.0828 

98153 

1.1490 

98172 

1.1883 

98402 

1.2333 

95297 

1.3117 

85250 

1.2866 

46917 

.9704 

49668 

1.0919 

64311 

1.1132 

83445 

1.1477 
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Figure Captions 

Fig. 1. Phase functions normalized to unity at four repre- 
sentative wavelengths in the near infrared for a cloud of ice 
spheres. The spheres follow the "cloud model" size distribution 
which has its maximum at diameter 8y. The inset on the right is 
an enlargement of the phase functions near the forward peak. 

Fig. 2. Spherical albedo in the near infrared of a plane 
parallel cloud of ice spheres with the "cloud model" size dis- 
tribution of particles. 

- 1 _ 0 

Fig. 3. Reflection function versus e = cos y for 4>-<{> 0 = 0 • 
The results are for the sharply peaked phase function typical of 
some clouds in the range 4000A-10 , 000A (the curve labeled X=.95 
in Fig. 1) . 

Fig. 4 . The same as Figure 3 but for <I>-<J> 0 = 60° 

Fig. 5. The same as Figure 3 but for <P~4> 0 - 180° 

Fig. 6. Reflection function versus azimuth angle 

<k— <k for & = cos - 1 y * 85°. The results are for the peaked 

T T 0 0 0 

phase function labeled X = .95y in Figure 1. The single scat- 
tering curves are for an optical thickness t = 1 but they are also 
approximately valid for thicker layers. 

A o 

Fig. 7. The same as Figure 6 but for e Q = 60 

o 

Fig. 8. The same as Figure 6 but for 0 Q = 30 
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